﻿using com.foxmail.wyyuan1991.Assignment.Solver.Data;
using System;
using System.Linq;

namespace com.foxmail.wyyuan1991.Assignment.Solver
{
    /// <summary>
    /// 对偶拉格朗日乘子法求解器
    /// </summary>
    public class SC_ALM_Solver
    {
        #region 属性
        /// <summary>
        /// 收敛阈值
        /// </summary>
        public double Epsilon { get; set; }//epsilon = 0.0000001m;
        /// <summary>
        /// 罚值
        /// </summary>
        public double Sigema { get; set; }//Sigema = 99m;
        /// <summary>
        /// 罚值/拉格朗日乘子 调整阈值
        /// </summary>
        public double Beta { get; set; }//beta = 0.001m;
        /// <summary>
        /// 罚值增加速度
        /// </summary>
        public double Alpha { get; set; }//alpha = 2m;
        /// <summary>
        /// 迭代上限
        /// </summary>
        public int IterationUpbound { get; set; }
        /// <summary>
        /// 最短路搜索算法
        /// </summary>
        public IRouteSearchAlgorithm SA { get; set; }
        /// <summary>
        /// 全有全无分配算法
        /// </summary>
        public AoNSolver AS { get; set; }
        #endregion

        #region Events
        //一次循环完成事件
        public event ALSIterationFinishedEventHandler AFEvent;
        private void RasieALSIterationFinishedEvent(int i,double p,double e)
        {
            if (AFEvent != null)
            {
                AFEvent(this, new ALSIterationFinishedEventArgs(i,p,e));
            }
        }
        #endregion

        #region 方法
        /// <summary>
        /// 对偶拉格朗日乘子法--Link Constraint
        /// </summary>
        /// <param name="net"></param>
        /// <param name="f"></param>
        /// <returns></returns>
        public IFlowPattern Solve(INetwork net, IDemand f)
        {
            CapacityLagrangeMultiplier lm = new CapacityLagrangeMultiplier();//这里的拉格朗日乘子是一维的，因为边约束的形式是link capacity.
            CapacityLagrangeMultiplier lm2 = new CapacityLagrangeMultiplier();
            //初始化拉格朗日乘子
            foreach (IArc a in net.ArcList.Where(i => i.Capacity > 0))
            {
                lm2[a]=lm[a] = a.Capacity;
            }
            //①全有全无分配
            FlowPattern fd,fd2 = AS.Solve(net, f) as FlowPattern;
            int iter = 1;
            //③罚值是否超限
            FWSolver s = new FWSolver()
                {
                    kesi = 0.000000001,//设置收敛阈值
                    SA = new DijstraRouteSearcher() { Network = net },
                    AS = this.AS,
                    LS = new DichotomyLineSearcher() { Kesi = 0.000000001 }
                };
            while (IsPenaltyLarge(net,lm,fd2) || iter == 1)
            {
                if (iter++ > IterationUpbound) break;
                fd = fd2; fd2 = null;
                //②内嵌F-W

                fd2 = s.Solve(net, f, delegate (IArc a, double load)
                {
                    return a.PerformaceFunction(load) + (a.Capacity > 0 ? (Math.Max(0, (lm2[a] - Sigema * (a.Capacity - load)))) : 0);
                }, fd) as FlowPattern;
                //触发单次迭代结束事件
                RasieALSIterationFinishedEvent(iter, CalculateTotalPenalty(net, lm2, fd2), CalculateTotalPenaltyEvolution(net,lm,lm2,fd,fd2));
                //④罚值改进幅度
                if (NeedtoEnlargeStep(net, lm, lm2, fd, fd2))
                {
                    Sigema = Alpha * Sigema;
                }
                UpdateLM(net, lm, lm2);
                //⑤更新乘子
                foreach (IArc a in net.ArcList.Where(i => i.Capacity > 0))
                {
                    lm2[a] = Math.Max(0, lm[a] - Sigema * (a.Capacity - fd2.getFlow(a)));
                }
            }
            return fd2;
        }
        //更新拉格朗日乘子
        private static void UpdateLM(INetwork net, CapacityLagrangeMultiplier lm, CapacityLagrangeMultiplier lm2)
        {
            foreach (IArc a in net.ArcList.Where(i => i.Capacity > 0))
            {
                lm[a] = lm2[a];
            }
        }
        //判断每条弧罚值是否超限
        private bool IsPenaltyLarge(INetwork net, CapacityLagrangeMultiplier lm, FlowPattern f)
        {
            foreach(IArc a in net.ArcList.Where(i=>i.Capacity>0))
            {
                if (calculatePenalty(a,f,lm) > Epsilon)
                {
                    return true;
                }
            }
            return false;
        }
        //计算总体罚值
        private double CalculateTotalPenalty(INetwork net, CapacityLagrangeMultiplier lm, FlowPattern f)
        {
            double sum = 0;
            foreach (IArc a in net.ArcList.Where(i => i.Capacity > 0))
            {
                sum += calculatePenalty(a, f, lm);
            }
            return sum;
        }
        //计算一条弧的罚值
        private double calculatePenalty(IArc a, FlowPattern f, CapacityLagrangeMultiplier lm)
        {
            double k =  Math.Abs(a.Capacity - f.getFlow(a) - (Math.Max(0, Sigema * (a.Capacity - f.getFlow(a)) - lm[a])/Sigema));
            return k;
        }
        //判断是否需要增大sigema
        private bool NeedtoEnlargeStep(INetwork net, CapacityLagrangeMultiplier lm, CapacityLagrangeMultiplier lm2,  FlowPattern f1, FlowPattern f2)
        {
            foreach (IArc a in net.ArcList.Where(i => i.Capacity > 0))
            {
                if (calculatePenalty(a, f1, lm)!=0&&calculatePenalty(a, f2, lm2)/calculatePenalty(a, f1, lm) > Beta)//防止除0错误，如果一条弧k-1次迭代中没有惩罚而k次有惩罚，只能说明方向不对，再增加sigama也没有用。
                {
                    return false;
                }
            }
            return true;
        }
        //计算总体罚值改进
        private double CalculateTotalPenaltyEvolution(INetwork net, CapacityLagrangeMultiplier lm, CapacityLagrangeMultiplier lm2, FlowPattern f1, FlowPattern f2)
        {
            if (CalculateTotalPenalty(net, lm, f1) != 0)
            {
                return CalculateTotalPenalty(net, lm2, f2) / CalculateTotalPenalty(net, lm, f1);
            }
            else
            {
                return -1;
            }
        }
        #endregion
    }
}
